1_IMR90_domains_23plex_analysis

Down stream analysis of the dataset ‘1_IMR90_domains_23plex’

setup

library(ggforce)
library(tidyverse)

dcolor <- ggsci::scale_color_d3('category20')
dfill <- ggsci::scale_fill_d3('category20')
ccolor <- viridis::scale_color_viridis()
style <- function(x,color,method='UMAP',axis=c(1,2),Theme=NULL,
                  axis.text=FALSE,axis.title=FALSE,coord_fix =TRUE,
                  palette=TRUE,title=NULL,
                  legend=TRUE,scale_color_log10=FALSE,
                  guide_pointSize=2.5,...) {
  axis <- paste0(method,axis)
  if(is.null(Theme)) Theme <- theme_bw
  g <- ggplot() + Theme()
  g <- g + geom_point(aes_(x=as.name(axis[1]),y=as.name(axis[2]),color=as.name(color)),x,...)
  
  if(coord_fix) {
    g <- g + coord_fixed()
  }
  if(!axis.text) {
    g <- g + theme(axis.text = element_blank(),axis.ticks = element_blank())
  }
  if(!axis.title) {
    g <- g + theme(axis.title = element_blank())
  }
  if(is.null(title)) {
    g <- g + ggtitle(method) 
  }else if(!is.na(title)) {
    g <- g + ggtitle(title) 
  }
  if(!legend) g <- g + theme(legend.position = 'none')
  
  if(scale_type(as.matrix(x[,color]))=='discrete' & !is.null(guide_pointSize)) {
    g <- g + guides(color = guide_legend(override.aes = list(size=guide_pointSize)))
  }
  
  if(palette) {
    if(scale_type(as.matrix(x[,color]))=='discrete') {
      g <- g + ggsci::scale_color_d3('category20')
    }else{
      if(scale_color_log10){
        g <- g + viridis::scale_color_viridis(trans='log10')
      }else{
        g <- g + viridis::scale_color_viridis()
      }
    }
  }
  return(g)
}
tib2df <- function(tib,RowNames=1) {
  rn <- pull(tib,all_of(RowNames))
  df <- tib %>%
    select(-RowNames) %>%
    as.data.frame()
  rownames(df) <- rn
  return(df)
}

load data

data <- tibble(
  path = list.files('1_IMR90_domains_23plex/Raw/',full.names = T),
  sampleName = paste(rep(c('G','S'),each=5),formatC(1:10,width = 2,flag = 0),sep = '-')
  ) %>%
  separate(sampleName,into=c('type','cellID'),remove = F) %>%
  mutate(metric = map(path,read_csv,show_col_types=F)) %>%
  select(-path) %>% unnest(metric) %>%
  unite(Domain,sampleName,domainID,sep = '-',remove = F) 

tmp <- data %>%
  gather(key=feature,value=exp,-(1:5)) %>%
  separate(feature,into=c("region","ab"),sep = " ") %>%
  mutate(region = factor(region,levels=c("In","Mid","Ex")),) %>%
  split(f=.$ab)
map(names(tmp),~{
  tmp[[.x]] %>%
    ggplot(aes(region,exp,color=cellID)) +
    theme_bw() + geom_violin(color="grey") + geom_jitter(size=.6) + 
    facet_wrap(~type) + ggtitle(.x) + dcolor +
    labs(y = 'Expression Level')
  }) %>%
  patchwork::wrap_plots(guides = "collect",ncol = 6)

Clustering

X <- as.matrix(tib2df(data[,-(2:5)]))
p <- prcomp(X)
screeplot(p,type='l',n=30)

total <- rowSums(data[,-(1:5)])
cor(p$x[,1:30],total) %>%
  plot(xlab = 'PC',ylab = 'CorCoef_vsTotalExp')

set.seed(10)
um <- uwot::umap(p$x[,2:10],metric = "cosine",pca = 10,n_neighbors = 10,min_dist = .3)
Meta <- data[,1:5] %>% 
  mutate(UMAP1=um[,1],UMAP2=um[,2])

style(Meta,"type",size=1,palette=F)

style(Meta,"cellID",size=1)

set.seed(1)
D <- as.matrix(dist(p$x[,2:10]))
k <- 10
adj <- apply(D,1,function(x) rank(x)%in%2:(k+1))
g <- igraph::graph_from_adjacency_matrix(adj,mode='undirected')
lei <- igraph::cluster_leiden(g,resolution_parameter = .02)

Meta <- Meta %>% mutate(cluster=as.factor(lei$membership))
style(Meta,"cluster")  + ggsci::scale_color_lancet()

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

tmp <- data %>% select(-sampleName,-type,-cellID,-domainID) %>%
  tib2df() %>% scale %>% t %>% scale %>% asinh()
idx <- rownames(tmp) %>% sub(".* ","",.)
mod <- model.matrix(~0+idx)
colnames(mod) <- sub('idx','',colnames(mod))
tmp2 <- (t(tmp)%*%mod)/3
Dists <- dist(t(tmp2))
h <- hclust(Dists, method = "average") 
dend <- as.dendrogram(h)
ggdend <- ggdendro::dendro_data(dend)
ggdend$segments$y <- ggdend$segments$y/3
ggdend$segments$yend <- ggdend$segments$yend/3
g <- ggplot(ggdendro::segment(ggdend))
g + geom_segment(aes(x=y,y=x,xend=yend,yend=xend)) + scale_x_reverse() + 
  geom_text(data = ggdend$label, aes(y, x, label = label), hjust = 0, angle = 0) + xlim(10,-8) + theme_void()

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

n <- ggdend$labels$label %>% length
idx <- ggdend$labels$label %>% rev %>% rep(each=3) %>% paste0(rep(c("In ","Mid ","Ex "),n),.) 
tmp <- tmp[idx,]

anno <- data %>% select(Domain,sampleName) %>% 
  separate(sampleName,into = c('type','cell')) %>% 
  tib2df()
anno$cluster <- as.factor(lei$membership)

anno.color <- list()
anno.color$type <- c('G' = "#F8766D", 'S' = "#00BFC4")
anno.color$cell <- ggsci::scale_color_d3()$palette(10)
anno.color$cluster <- ggsci::scale_color_lancet()$palette(4)
names(anno.color$cell) <- formatC(1:10,width = 2,flag = '0')
names(anno.color$cluster) <- 1:4

pheatmap::pheatmap(tmp[,lei$membership == 1],
                   cluster_rows = F,
                   annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
                   color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)

pheatmap::pheatmap(tmp[,lei$membership == 2],
                   cluster_rows = F,
                   annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
                   color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)

pheatmap::pheatmap(tmp[,lei$membership == 3],
                   cluster_rows = F,
                   annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
                   color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)

pheatmap::pheatmap(tmp[,lei$membership == 4],
                   cluster_rows = F,
                   annotation_col = anno,annotation_colors = anno.color,show_colnames = F,
                   color = viridis::viridis(32),border_color = NA,cellheight = 10,cellwidth = 2.5)

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] SeuratObject_4.1.3 Seurat_4.3.0.1     ggforce_0.4.1      lubridate_1.9.2    forcats_1.0.0      stringr_1.5.0      dplyr_1.1.2        purrr_1.0.1        readr_2.1.4       
[10] tidyr_1.3.0        tibble_3.2.1       ggplot2_3.4.2      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] utf8_1.2.3             spatstat.explore_3.2-1 reticulate_1.30        tidyselect_1.2.0       RSQLite_2.3.1          AnnotationDbi_1.60.2   htmlwidgets_1.6.2     
  [8] grid_4.2.2             BiocParallel_1.32.6    Rtsne_0.16             munsell_0.5.0          codetools_0.2-19       ragg_1.2.5             ica_1.0-3             
 [15] future_1.33.0          miniUI_0.1.1.1         withr_2.5.0            spatstat.random_3.1-5  colorspace_2.1-0       progressr_0.13.0       Biobase_2.58.0        
 [22] knitr_1.43             rstudioapi_0.15.0      stats4_4.2.2           ROCR_1.0-11            tensor_1.5             listenv_0.9.0          labeling_0.4.2        
 [29] GenomeInfoDbData_1.2.9 polyclip_1.10-4        bit64_4.0.5            farver_2.1.1           pheatmap_1.0.12        parallelly_1.36.0      vctrs_0.6.3           
 [36] generics_0.1.3         xfun_0.39              timechange_0.2.0       R6_2.5.1               phateR_1.0.7           GenomeInfoDb_1.34.9    bitops_1.0-7          
 [43] spatstat.utils_3.0-3   cachem_1.0.8           promises_1.2.0.1       scales_1.2.1           vroom_1.6.3            gtable_0.3.3           globals_0.16.2        
 [50] goftest_1.2-3          rlang_1.1.1            systemfonts_1.0.4      RcppRoll_0.3.0         pamr_1.56.1            splines_4.2.2          lazyeval_0.2.2        
 [57] spatstat.geom_3.2-4    yaml_2.3.7             reshape2_1.4.4         abind_1.4-5            stylo_0.7.4            httpuv_1.6.11          tools_4.2.2           
 [64] tcltk_4.2.2            ellipsis_0.3.2         jquerylib_0.1.4        RColorBrewer_1.1-3     ggdendro_0.1.23        proxy_0.4-27           BiocGenerics_0.44.0   
 [71] ggridges_0.5.4         Rcpp_1.0.11            plyr_1.8.8             zlibbioc_1.44.0        RCurl_1.98-1.12        deldir_1.0-9           pbapply_1.7-2         
 [78] viridis_0.6.4          cowplot_1.1.1          S4Vectors_0.36.2       zoo_1.8-12             ggrepel_0.9.3          cluster_2.1.4          magrittr_2.0.3        
 [85] data.table_1.14.8      scattermore_1.2        lmtest_0.9-40          RANN_2.6.1             fitdistrplus_1.1-11    Signac_1.10.0          matrixStats_1.0.0     
 [92] hms_1.1.3              patchwork_1.1.2        mime_0.12              evaluate_0.21          xtable_1.8-4           IRanges_2.32.0         gridExtra_2.3         
 [99] compiler_4.2.2         KernSmooth_2.23-22     crayon_1.5.2           htmltools_0.5.5        later_1.3.1            tzdb_0.4.0             DBI_1.1.3             
[106] tweenr_2.0.2           formatR_1.14           MASS_7.3-60            Matrix_1.5-4.1         cli_3.6.1              parallel_4.2.2         igraph_1.5.0.1        
[113] GenomicRanges_1.50.2   pkgconfig_2.0.3        sp_2.0-0               plotly_4.10.2          spatstat.sparse_3.0-2  bslib_0.5.0            XVector_0.38.0        
[120] digest_0.6.33          sctransform_0.3.5      RcppAnnoy_0.0.21       tsne_0.1-3.1           spatstat.data_3.0-1    Biostrings_2.66.0      rmarkdown_2.23        
[127] leiden_0.4.3           fastmatch_1.1-3        uwot_0.1.16            shiny_1.7.4.1          Rsamtools_2.14.0       lifecycle_1.0.3        nlme_3.1-162          
[134] jsonlite_1.8.7         viridisLite_0.4.2      fansi_1.0.4            pillar_1.9.0           ggsci_3.0.0            lattice_0.21-8         KEGGREST_1.38.0       
[141] fastmap_1.1.1          httr_1.4.6             survival_3.5-5         glue_1.6.2             png_0.1-8              bit_4.0.5              tcltk2_1.2-11         
[148] class_7.3-22           stringi_1.7.12         sass_0.4.7             blob_1.2.4             textshaping_0.3.6      memoise_2.0.1          irlba_2.3.5.1         
[155] e1071_1.7-13           future.apply_1.11.0    ape_5.7-1